package fem2.material;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import math2.LinAlgUtilities;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.NotImplementedError;
import fem2.element.StructuralElement;
import fem2.enu.IntegrationType;
import fem2.enu.PropertyName;
import fem2.enu.State;

public abstract class ElasticMaterialModel extends MaterialModel {

	protected State state;
	protected double kappa;
	protected double lambda;
	protected double mu;

	public ElasticMaterialModel(double E, double nu, double thickness, double rho, State ss) {
		super(rho, thickness, ss);

		kappa = E / (3 * (1 - 2 * nu));
		mu = E / (2 * (1 + nu));
		lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
		this.state = ss;
	}

	@Override
	public double getMaterialProperty(PropertyName name) {
		switch (name) {
		case SHEAR_MODULUS:
			return mu;
		case BULK_MODULUS:
			return kappa;
		case POISSON_RATIO:
			return lambda / (2 * (lambda + mu));
		case YOUNG_MODULUS:
			double nu = getMaterialProperty(PropertyName.POISSON_RATIO);
			return mu * 2 * (1 + nu);
		default:
			throw new NotImplementedError();
		}
	}

	/**
	 * compute strain energy term given strain at material point
	 * 
	 * @param e
	 * @return
	 */
	public abstract double computeEnergyAt(IVector e);

	/**
	 * compute strain energy term given stress & strain at material point
	 * 
	 * @param e
	 * @return
	 */
	public abstract double computeEnergyAt(IVector e, IVector s);

	/**
	 * create elasticity matrix for plane stress state
	 * 
	 * @param e
	 *            elastic modulus
	 * @param nu
	 *            Poisson ratio
	 * @return De matrix
	 */
	public static IMatrix createPlaneStressElasticityMatrix(double e, double nu) {
		IMatrix c = Array1DMatrix.createRowMajor(3, 3);
		double k = e / (1 - Math.pow(nu, 2));
		c.set(0, 0, k * 1);
		c.set(1, 0, k * nu);
		c.set(0, 1, k * nu);
		c.set(1, 1, k * 1);
		c.set(2, 2, k * (1 - nu) / 2);
		return c;
	}

	/**
	 * create elasticity matrix for plane strain state
	 * 
	 * @param e
	 *            elastic modulus
	 * @param nu
	 *            Poisson ratio
	 * @return De matrix
	 */
	public static IMatrix createPlaneStrainElasticityMatrix(double e, double nu) {
		IMatrix c = Array1DMatrix.createRowMajor(3, 3);
		double k = e / ((1 + nu) * (1 - 2 * nu));
		c.set(0, 0, k * (1 - nu));
		c.set(1, 0, k * nu);
		c.set(0, 1, k * nu);
		c.set(1, 1, k * (1 - nu));
		c.set(2, 2, k * (1 - 2 * nu) / 2);
		return c;
	}

	/**
	 * create three dimensional elasticity matrix
	 * 
	 * @param e
	 *            elastic modulus
	 * @param nu
	 *            Poisson ratio
	 * @return De matrix
	 */
	public static IMatrix create3dElasticityMatrix(double e, double nu) {
		IMatrix c = Array1DMatrix.createRowMajor(6, 6);
		double k = e / ((1 + nu) * (1 - 2 * nu));
		double c1 = k * (1 - nu);
		double c2 = k * nu;
		double c3 = e / (2 * (1 + nu));

		c.set(0, 0, c1);
		c.set(0, 1, c2);
		c.set(0, 2, c2);

		c.set(1, 0, c2);
		c.set(1, 1, c1);
		c.set(1, 2, c2);

		c.set(2, 0, c2);
		c.set(2, 1, c2);
		c.set(2, 2, c1);

		c.set(3, 3, c3);
		c.set(4, 4, c3);
		c.set(5, 5, c3);
		return c;
	}

	/**
	 * compute strain energy of this element
	 * 
	 * @return strain energy
	 */
	public double computeStrainEnergy(StructuralElement se) {
		MeshPart mp = se.getMeshPart();
		double energy = 0.0;
		double th = getSectionProperty();
		int ng = se.getNumIntegrationPoints();

		for (int k = 0; k < ng; k++) {
			double w = se.getIntegrationWeight(k);
			double[] xi = se.getIntegrationPoint(k);

			double[][] jac = mp.jacXAt(xi);
			double detJ = LinAlgUtilities.determinantNxN(jac);
			double coeff = w * detJ * th;

			IVector e = se.computeStrain(xi);
			IVector s = se.getStress(k);
			energy += computeEnergyAt(e, s) * coeff;
		}
		if (se.getIntegrationType() == IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			energy /= 2;
		}
		return energy;
	}
}
